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ABSTRACT 

We present results from the first high-resolution hydrodynamical simulations of non- 
Gaussian cosmological models. We focus on the statistical properties of the transmitted 
Lyman-a flux in the high redshift intergalactic medium. Imprints of non-Gaussianity 
are present and are larger at high redshifts. Differences larger than 20% at z > 3 in 
the flux probability distribution function for high transmissivity regions (voids) are 
expected for values of the non linearity parameter /nl = ±100 when compared to a 
standard ACDM cosmology with /nl = 0. We investigate also the one-dimensional flux 
bispectrum: at the largest scales (corresponding to tens of Mpc) we expect deviations 
in the flux bispectrum up to 20% at z ~ 4 (for /nl = ±100), significantly larger 
than deviations of ~ 3% in the flux power spectrum. We briefly discuss possible 
systematic errors that can contaminate the signal. Although challenging, a detection 
of non-Gaussianities in the interesting regime of scales and redshifts probed by the 
Lyman-a forest, could be possible with future data sets. 
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1 INTRODUCTION 

According to the standard gravitational instability pic- 
ture present-day cosmic structures have evolved from tiny 
initial fluctuations in the mass density field that obey 
Gaussian statistics. However, departures from Gaussian- 
ity inevitably arise at some level during the inflationary 
epoch. The various mechanisms that produce primordial 
non-Gaussiani ty during inflation have been thoroughly in- 
vestigated bv lBartolo et al. (2004) (and references therein). 
A convenient way of modeling non-Gaussianity is to include 
quadratic correction in the Bardeen's gauge-invariant poten- 
tial $: 

<£> = $ L + /nl ($l - ($£)) , (1) 

where $l represents a Gaussian random field and the di- 
mensionless parameter /nl quantifies the amplitude of the 
corrections to the curvature perturbations. The above defi- 
nition in which the term — /nl (<&l) i s small guarantees that 
< $ > = < <&l >= 0. Although the quadratic model quan- 
tifies the level of primordial non-Gaussianity predicted by 



a large number of scenarios for the generation of the ini- 
tial seeds for structure formation (including standard single- 
field and multi-field inflation, the curvaton and the inho- 
mogeneous reheating scenarios), one should keep in mind 
that there are different ways for a density field to be non- 
Gaussian (NG) and that different observational tests capa- 
ble of going beyond second order statistics should be used 
to fully characterize the nature of non-Gaussianity. 

To date, the strongest observational constraint for NG 
models are provided by the recent a nalysis of the WMAP 
5-year temperature fluctuation maps (|Komatsu et al.ll2008l ) 
according to which —9 < /nl < 111 at the 95 % confidence 
level in the local model. The large scale structure (LSS) 
provides alternative observational constraints which are, in 
principle, more stringent than the cosmic microwave back- 
ground (CMB) since they carry information on the 3D pri- 
mordial fluctuation fields, rather than on a 2D temperature 
map. Moreover, if the level of primordial non-Gaussianity 
depends on scale then CMB and LSS provide independent 
constraints since they probe different scales. For this reason, 
the WMAP 5-year limits on /nl need not to apply on the 
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smaller scales probed by the LSS and the NG models that 
we consider in this work, which have |/nl| as large as 200, 
are thus not in conflict with the CMB on the scales which 
are relevant for our analysis. 

A very promising way to constrain departures 
from Gaussianity is to measure the various proper- 
ties of massive virializ e d structures like t heir abundance 
dMatarrese et all |2000| ; IVerde et all l200ll; ILo Verde et all 
20081). clustering and their biasing jGrinstein fc Wisd ll986l; 
Matarrese et al. 1986; Matarrese fc Verde 2008; Dalai et al 



2008 ; ICarbone et al-feoosl ; ISeliakl 120081 : I Matarrese fc Verde 
20081 ). Indeed, the best constr aints on non-Gaussi anity from 



the LSS have been obtained bv lSlosar et al.l l|2008h including 
the observed scale-dependent bias of the spectroscopic sam- 
ple SDSS luminous red galaxies and the photometric quasar 
sample. The resulting limits of —29 < /nl < 70 (95 % confi- 
dence level) are remarkably close to tho se obta i ned fr om the 
CMB analysis alone and, according to ISeliakl l|2008h . could 
be further improved by looking for scale dependency in the 
relative biasing of two different population of objects. Alter- 
nativ ely, one can consi der the topology of the mass density 
field (jMatsubara 200jl), and higher-order c lustering statis- 
tics like the bispectrum (jHikage et al.l l2006). The ability of 
these techniques to detect the imprint of the primordial 
non-Gaussianity on the LSS has been tested with N-body 



experiments dMessina et al.l 


1990; 


1 

Moscardini ct al. 


11991 


Weinberg & Cole 1992; Mathis et al. 2004: Kang et al 


2007; 


Grossi ct al. 2007; Dalai ct al. 2008; Hikagc et al. 2008). N- 



body simulations are of paramount importance in the study 
of NG models, since one needs to disentangle primordial 
non-Gaussianity from late non-Gaussianity induced by the 
non-linear growth of density perturbations that can only be 
properly accounted for by numerical experiments. 

Recently, iGrossi et al. I (|2007l .l 2008) have carried out cos- 
mological N-body simulations of NG models to study the 
evolution of the probability distribution function (PDF) of 
the density fluctuations. They found that the imprint of 
primordial non-Gaussianity, which is evident in the nega- 
tive tail of the PDF at high redshifts, is preserved through- 
out the subsequent evolution and out to the present epoch. 
This result suggests that void statistics may be a promis- 
ing effective tool for detect in g primordial non-G aussianity 
i|Kamionkowski et all (|200Sf ): ISong fc Leel ((2008)) and that 
can be applied to different types of observations over a large 
range of cosmic epochs. Taking advantage of the recent the- 
oretic al efforts for s tandardizing the appropriate statistical 
tools (Colbcrg 2008) one could apply void-finding algorithms 
to quantify the properties of the underdense regions ob- 
served in the spatial distribution of galaxies. Unfortunately, 
current galaxy redshift surveys are probably too small for 
void-based statistics to appreciate deviations from the Gaus- 
sian case at the level required. The situation will change in a 
not too distant future, when next generation all-sky surveys 
like ADEPT or EUCLID will allow to measure the position 
of ~ 5 x 10 r galaxies over a large range of redshifts out to 
z = 2. Alternatively, one can analyze high-resolution spec- 
tra of distant quasars to characterize the properties of the 
unde rlyi ng mass density field at z > 3 fe.g. lViel etHI (|2003l 
l2004l lbl); iLesgourgues et al] l|2007l )). In particular, since we 
expect that underdense regions are characterized by a low 
neutral hydrogen (HI) abundance, one can infer the presence 
of voids and quantify their statistical properties from voids 



in the transmitted flux, defined as the connected regions in 
the spectral flux distribution above the mean flux level. The 
connection between voids and spectral regions characterized 
by negligi b le HI absorption has been recently studied by 
IViel et all (|20Qgft using hydrodynamical simulations where 
a link at z ~ 2 between the flux and matter properties is 
provided. 

In this work we perform, for the first time, high- 
resolution hydrodynamical simulations of NG models to 
check whether one can use the intergalactic medium (see 
iMeiksinl (|2007T ) for a recent review) to detect non-Gaussian 
features in the Lyman-a flux statistics like the PDF, flux 
power and the bispectrum. The layout of the paper is as fol- 
lows. In Section 2 we describe the hydrodynamical simula- 
tions and we show an example of simulated Lyman-a quasar 
(QSO) spectrum. In Section 3 we present the results of the 
various flux statistics. In Section 4 we address the role of 
systematic and statistical errors that could contaminate the 
NG signal. We conclude in Section 5. 



2 NON-GAUSSIAN HYDRODYNAMICAL 
SIMULATIONS 

We rely on simulations run with the parallel hydrodynam- 
ical (TreeSPH) code GADGE T-2 based on th e conservative 
'entropy-formulation' of SPH (Springel 2005). They consist 
of a cosmological volume with periodic boundary conditions 
filled with an equal number of dark matter and gas parti- 
cles. Radiative cooling and heating processes were followed 
for a primordial mix of hydrogen and helium. We assumed 
a mean Ultraviole t Back ground similar to that propesed by 
lHaardt fc~ Madau (1996) produced by quasars and galaxies 
as given by with helium heating rates multiplied by a fac- 
tor 3.3 in order to better fit observatio nal constraints on th e 
temperature e volut ion of the IGM (e.g. lSchave et al.l (2000); 
iRicotti et all (|200d lV This background gives naturally a hy- 
drogen ionizat ion rate T-i2 ~ 1 at the redshifts of inter- 
est he re (e.g. iBolton et alj (|2005l ); iFaucher-Giguere et al.l 
(2008)). The star formation criterion is a very simple one 
that converts in collisionless stars all the gas particles whose 
temperature falls below 10 5 K and whose density contrast is 
larger than 1000 (it has been shown that the star formation 
criterion has a negligib le impact on flu x statistics). More 
details can be found in j|Viel et alj|2003 ). 

The cosmological reference model corresponds to a 
'fiducial' ACDM Universe with parameters, at z = 0, f2 m = 
0.26, fl A = 0.74, fl b = 0.0463, n s = 0.95, an d Hp = 72 
km s" 1 Mpc" 1 and erg = 0.85 (the B2 series of I Viel et ail 
(2004)). We have used 2 x 384 3 dark matter and gas particles 
in a 60 comoving Mpc box for the flux power and bispec- 
trum, to better sample the large scales. For the flux proba- 
bility distribution function we relied instead on 2 x 256 3 dark 
matter and gas particles in a 20 ft -1 comoving Mpc, since 
below and around z — 3 this seems to be the appropriate 
resolution the get numerical convergence. The gravitational 
softening was set to 2.5 and 5 h -1 kpc in comoving units 
for all particles for the 20 and 60 comoving Mpc/ft boxes, 
respectively. The mass per gas particle is 6.12 x 1O 6 M0 /h for 
the small boxes and 4.9 x 10 7 Mq jh for the large boxes, while 
the high resolution run for the small box has a mass per gas 
particle of 1.8 x 1O 6 M //i (this refers to a (20,384) simu- 
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Figure 1. Gas (projected) density distribution in the hydrodynamical simulations at z = 3 in the Gaussian case (middle panel) and 
residuals in the gas distribution for non-Gaussian models with /nl = —200 (left panel) and /nl = +200 (right panel). The simulated 
box has a linear size of 20 comoving Mpc//i and the thickness of the slice shown is 1.5 comoving Mpc//i. The cosmic web that gives 
rise to Lyman-a absorption is visible and the tendency for regions below the mean density to be underdense (overdense) in models with 
negative (positive) values of /nl is clear. 



lation that was performed in order to check for numerical 
convergence of the flux PDF). In the following, the different 
simulations will be indicated by two numbers, (Ni, N2): Ni 
is the size of the box in comoving Mpc/ft and N2 is the cubic 
root of the total number of gas particles in the simulation. 
NG are produced in th e initial cond i tions at z = 99 using 
the same method as in iGrossi et all (|2007f ) that we briefly 
summarize here. Initial NG conditions are generated with- 
out modifying the linear matter power spectrum using the 
Zel'dovich approximation: a Gaussian gravitational poten- 
tial is generated in Fourier space from a power-law power 
spectrum of the form P(k) oc k~ 3 and inverse Fourier trans- 
formed in real space to produce </>l- The final <E> is obtained 
using eq. (1). Finally, back in Fourier space, we modulate 
the power-law spectrum using the transfer functions of the 
ACDM model. 

We also run some other simulations at higher resolu- 
tions to check for numerical convergence. In particular we 
have performed a (20, 384) simulation run to analyse the flux 
PDF. For the (20, 256) models the flux probability distribu- 
tion fun£tion_Jias nume rically converged only below z = 3 
(see lBolton et all l|2008h ). However, since our results will al- 
ways be quoted in comparison with the /nl = case (i.e. as 
a ratio of two different quantities) we expect the resolution 
errors to be unimportant (i.e. we assume the same resolution 
corrections should be applied to all the models, even though 
this assumption should be explicitly checked). 

A projected density slice of the gas (IGM) distribution 
for the (20,256) simulation of thickness 1.5 comoving Mpc/h 
is shown in Figure [T] We focus on this simulation because at 
z — 3 the flux probability distribution function has numeri- 
cally converged. In the middle panel we plot the gas density 
in the Gaussian case, while residuals in the two models with 
/nl = —200 and 5ml = +200 are shown in the left and right 
panel, respectively. On average regions of the cosmic web 
below the mean density appear to be ~ 10% less (more) 
dense in the negative (positive) /nl case. This trend is ap- 



parent not only near the centre of these regions but also in 
the matter surrounding them (see for example the void at 
(x — 17; y — 8) comoving Mpc/h). The same qualitative be- 
havior can be observed in the distribut i on of t he dark matter 
particles (see Figure 2 of lGrossi etafl (|2008h V 

In the NG models considered here the growth of struc- 
tures in terms of density PDF is different. As discussed 
in Grossi et al. (2008), the maps of residuals in the non- 
Gaussian cases reflect the differences in the primordial PDF 
of the mass overdensity. As shown in Figures 1 and 5 of 
Grossi et al. (2008) the mass PDF is skewed towards positive 
(negative) overdensities in the non-Gaussian models with 
positive (negative) /nl values, compared to the Gaussian 
case. As a consequence, since the gas traces well the under- 
lying mass distribution at these redshifts, voids look emptier 
in the /nl = —200 case (map on the left) while denser en- 
vironments like filaments and knots look more prominent in 
the /nl = +200 case (map on the right) with respect to the 
Gaussian case. These differences in the tails of the density 
PDF also impact on the filaments at around the mean den- 
sity that surround the voids. In fact the size of the voids is 
slightly different in the negative and positive NG models: 
for negative /nl values the emptier voids grow in size faster 
than for the Gaussian case and even faster than for positive 
/nl values, displacing the filaments at around the mean den- 
sities at different positions in the three cases and giving rise 
to the filamentary pattern of residuals of the panels. 

To perform our analysis we have extracted several mock 
QSO absorption spectra from the simulation box. All spectra 
are drawn in redshift space taking into account the effect 
of the IGM peculiar velocities along the line-of-sight v peCi u. 
Basically, the simulated flux at the redshift-space coordinate 
u (in km/s) is F(u) = exp[— t(u)] with: 

/oo 
dxnm(x) Q[u-x- v%Z™\(x), b(x)] dx , (2) 

where ao, a = 4.45 x 10 -18 cm 2 is the hydrogen Lyman- 
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Figure 2. Simulated noiseless Lyman-ct flux (only regions at high transmissivity are shown) at z = 3 in the upper panel. Models with 
positive and negative values of /nl bracket the Gaussian case with /jjl = 0. Although difficult to see on a pixel-by-pixel basis, differences 
among the models can be appreciated in a statistical sense when many spectra are considered (see text). The bottom panel shows the 
ratio of non-Gaussian models to the Gaussian one for the one-dimensional gas density along the line-of-sight (in real space). On average 
differences of the order 10% are present. 



q cross-section, H(z) is the Hubble constant at red- 
shift 2, x is the real-space coordinate (in km s -1 ), b = 
(2fcflT/mc 2 ) 1 ' 2 is the velocity dispersion in units of c, 
9 = (0F6)- 1 exp[-( U - y - vl™(y)) 2 /b 2 ] is the Gaus- 
sian profile that well approximates the Voigt profile in the 
regime considered here. The neutral hydrogen density in 
real-space, that enters the equation above, could be related 
to t he underlying gas de n sity by t he fo llowing expression 
(e.g. iHui fc GnedhTll 19971 ); ISchavd l|200ll )): 



nm(x, z) 



10- 5 n IGM (*) 



■ T(x,z) ' 
v 10 4 K , 



(i^) 3 (l + 5lGM(x,2)) 2 



(3) 



with T-12 is the hydrogen photoionization rate in units of 
s _1 , T is the IGM temperature, and n\Gm{z) is the mean 
IGM density at that redshift. However, this equation is not 
explicitly used since the neutral hydrogen fraction is com- 
puted self-consistently for each gas particles during the sim- 
ulation run. The integral in eq. @ to obtain the Lyman- 
a optical depth along each simulated line-of-sight is thus per- 
formed using the relevant hydrodynamical quantities from 
the numerical simulations: 5igm, T, v pec , nm- More details 
on how to extract a mock QSO spectrum from an hydrody- 
n amical simulation usi ng the SPH formalism can be found 
in lTheuns et all (|l998l ). 

An example of line-of-sight is shown in the top panel of 
Figure [2] while the bottom panel shows the ratio of the gas 
density along the line-of-sight of NG and Gaussian mod- 
els (in real space) In the following we will focus on high 



transmissivity in which the transmitted flux is close to unity 
(upper panel). Three QSO spectra are shown with different 
line styles and correspond to the Gaussian case (dashed, red 
line), and to /nl = ±200 (solid black and solid blue, re- 
spectively). The transmitted flux (no noise is added in this 
case) is almost identical for the two NG models in magni- 
tude but not in sign, as expected. One can better appreci- 
ate the differences among the models by looking at the gas 
density (bottom panel). On average differences are of the 
order 10%, even if in some cases they can rise above 30- 
40%. The fact that the corresponding variations in the flux 
are comparatively smaller (usually less than few percent) is 
somehow expected, since differences in the gas density are 
exponentially suppressed by the non-linear transformation 
between flux and matter (and by other non-linear effects as 
well). However, despite their small amplitude, the differences 
in the transmitted flux are large enough to be appreciated 
through appropriate statistical analyses of many indepen- 
dent lines-of-sight, as we will see in the following sections. 
Global statistics will be usually shown for samples of 1000 
lines-of-sight extracted along random directions within the 
simulated volume. 
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Figure 3. Ratio between the simulated flux probability distribution function of four different (20, 256) models with /nl = 
—200, —100, +100, +200 and /nl = 0, represented by the blue dashed, blue continuous, black continuous and black dashed lines, 
respectively. Results are shown at z = 3,4,5.5 in the left, middle and right panels. (Note the different scale in the right panel). 



3 RESULTS 



3.1 The Lyman-a flux probability distribution 
function 

In Figure [3] we show results for the flux probability distribu- 
tion function at z — 3, 4 and 5.5 using the (20,256) simula- 
tions, in the left, middle and right panel, respectively. The 
mock QSO spectra have been normalized to reproduce the 
same (observed) mean flux level. The scaling factor is usu- 
ally different by less than 2% to the standard non-Gaussian 
case (more precisely the differences are below 1% at z~3, 
around 1% at z=4.0 and around 2% at z=5.5 between the 
Gaussian and the /nl = ±200 cases). Differences between 
the non-Gaussian case and the Gaussian one are apprecia- 
ble only in regions of high transmissivity (flux ~ 1), that are 
typically associated to connected regions below mean den- 
sity (voids) in the matter distribution. At z — 3 (left panel) 
the differences can be of the order of 20% (40%) for models 
with /nl = ±100(±200). Models with negative (positive) 
values of /nl produce less (more) absorption. This reflects 
the fact that voids in models with negative /nl are emptier 
of neutral hydrogen than in the Gaussian case. The opposite 
holds true for models with positive /nl- T his is analogous to 
the effect discussed bv lGrossi et al.l l|2008l ) on the dark mat- 
ter density field and characterized in terms of the probability 
distribution function of density fluctuations. In that case for 
negative values of /nl the low density tail of the dark matter 
density PDF is more prominent. In our case what is more 
prominent is the high flux tail of the Lyman-a flux PDF. 
The amplitude of the effect increases with the redshift. At 
z — 4 (middle panel) differences w.r.t. the Gaussian case 
are as large as 30%-60% (for /nl = ±100 and ±200, re- 
spectively) and at z = 5.5 (right panel) the differences are 
of the order of ~ 100 % - 150% (for /nl = ±100 and ±200, 
respectively). Note that in the latter (z = 5.5) case we have 
used a different scale for the y— axis. 

From an observational viewpoint, it should be noted 



that the Lyman-a flux PDF has been measured with great 
accuracy using high-resolution spectra taking into account 
the metal contaminati ons and continuu m fitting errors at 
z = 2.07,2.54,2.94 by iKim et all (|2007l ). On the contrary, 
continuum fitting errors and the metal contaminations are 
somewhat harder to estimate in th e measuremen t s at h igher 
redshifts (z = 4.5 and z = 5.5) by iBecker etHI l|2007h . We 
will come back to this point in Section 4. 



3.2 The Lyman-a flux void distribution function 

A different, although not completely unrelated statistics 
is represented by the probability distribution function of 
the voids of given comoving size R. Searching for voids in 
the Lyman-a forest of observed QSO spectra has a lon g 



1987) 



1988 



dating histor y (see for example Carswell fc Reea 
Crott j (Il987l); iDuncan et alJ <ll989l): lOstriker et all v 
Dobrzvcki fc Bechtoldl l| 19911 ): iRauch et alJ l|l992l )) but in 
this paper we focus on the impact of non-Gaussianities on 
their statistical properties. 

We define flux voids as in IViel et all <|2008l ): connected 
one-dimensional regions along the QSO spectrum whose 
transmitted flux is above the mean flux level at that redshift. 
In Fig. |4]we show the ratio between the probability distribu- 
tion functions of the NG and Gaussian cases at z — 3.0. For 
this plot we have used the (60,384) that have the largest box 
size. Although the size of the largest voids (R~ 20 comoving 
Mpc/ft) is comparable to that of the box the corresponding 
differences in the probability distribution are rather mild 
and of the order of 10-15%. 

The differences are as expected for voids of sizes larger 
than 20 comving Mpc/ft (while for smaller ones the differ- 
ences are negligible): negative values of /nl result in voids 
that are emptier compared to the standard Gaussian case 
and thereby the typical sizes could be larger; while the op- 
posite trend can be seen for the positive values of /nl- The 
effect, even for /nl = ±200, is however somewhat smaller 
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Figure 4. Ratio of the fraction of voids in the flux distribu- 
tion between NG models and the corresponding Gaussian one at 
z = 3 as a function of void size in comoving Mpc/h. The models 
with / NL = -200, -100, +100, +200 are represented by the blue 
dashed, blue continuous, black continuous and black dashed lines, 
respectively. 



than the effects that can be induced by changing other cos- 
m ological or astroph ysical parameters (see the relevant plots 
in IViel et al.l (|2008l )). Furthermore, the uncertainty in the 
mean flux level at z = 3, which enters the definition of a 
void in the flux, produces an effect that is still larger than 
the NG signal sought. 

3.3 The Lyman-a flux power spectrum 

Primordial non-Gaussianity affects the evolution of density 
perturbation, particularly at the epochs and scales in which 
they enter the non-linear regime. Deviations from the Gaus- 
sian case are larger at high redshift, since at late times 
nonlinear dynamical effects become dominant. However, the 
contribution non-Gaussianity implied by /nl = ±100 is al- 
ways within a few percent of the total gravitational potential 
and should not appreciably affect the linear matter power 
spectrum. Consequently, we also expect the effect on the ID 
flux power spectrum to be small. 

To quantify the effect we have plotted in Figure [5] the 
ID flux power spectrum for the Gaussian and non-Gaussian 
cases at z = 3, 4, and 5.5, in the left, middle and right panels 
respectively. Even in this case the QSO spectra have been 
normalized to reproduce the same mean flux. Differences to 
the Gaussian case are of the order of 2%, 3% and 5% at 
the redshifts considered here and manifest themselves as an 
overall plateau with slightly more power at the largest scales 
(a factor two larger than at the smallest scales probed). As 
expected, the effect of primordial non-Gaussianity on the 
flux power spectrum is small and the effect decreases with 
time. 

In principle this effect on the flux power is degener- 
ate only with a change in the mean flux level (see for ex- 



ample Figure 3 o f I Viel fc Haehneltl (120061 ) or Figure 13 of 
iMcDonald et~aH (|2005l l'l: this means that other changes in 
cosmological parameters and/or astrophysics produce a dif- 
ferent fc-dependent change in the flux power than the one 
produced by non-Gaussianities. However, the magnitude of 
this effect is quite small and probably not detectable with 
present data sets. 



3.4 The Lyman-a flux bispectrum 

Unlike the power spectrum, the bispectrum on large scales 
is sensitive to the statistical properties of primordial fluctu- 
ation s like a primordial non - Gaussianity l|Fr^lll994l ; IVerdel 
120021 : ISefusatti fc Komatsul [20071 ) . Therefore, the ID flux 
bispectrum looks like a very promising statistics to search 
for non-Gaussianities in the IGM. The Lyman-a flux bis- 
pectrum has been calculat ed for the first tim e using high- 
resolution QSO spectra bv lViel et al] (|2004bh . Here, we use 
the same definition i.e. the real part of the three point func- 
tion in fc— space ,Df = Re(<5p(fci) 5f(&2) ^(^3)), for closed 
triangles fci + ki + kz = 0. 5p(fc) is the Fourier transform of 
5F. Df is related to the bispectrum of the flux Bp(ki, k2,kz) 



(D F ) =2nB F {k 1 ,k 2 ,k 3 )5 D (k 1 + k2 + k 3 ) 



(4) 



5 D (k) is the one-dimensional Dirac delta function and (■} 
indicates the ensemble average. Since we compute the one- 
dimensional bispectrum our triangles are degenerate and we 
choose two configurations: i) the flattened configurations for 
which fci = k-i and kz = — 2fci; ii) the squeezed configura- 
tion for which fci = fc — fc m i n , /c2 = — fc — fcmin and kz = 2 fcmin, 
with fc m i n — 2n/L (L the linear size of the box in km/s). In 
the following we will always show the flux bispec t rum as a 
function of the wavenumber k = k%. In I Viel et al] (|2004bh a 
numerical calculation of the flux bispectrum was compared 
to analytical estimates obtained through an expansion at 
seco nd order of the fluctuatin g Gunn-Peterson approxima- 
tion (|Gunn fc Peterson! (|l965T )V while the overall amplitude 
of the bispectrum was not matched by the theory, the shape, 
at least at large scales, was well reproduced. However, the 
theoretical expression for the flux bispectrum contained only 
the gravitational terms. Here we extend this work by com- 
puting the flux bispectrum for NG Gaussian models using 
the numerical hydrodynamical simulations performed. 

In Figure |6] we plot our findings in terms of ratios 
between the Gaussian and non-Gaussian models in the 
squeezed (top panels) and flattened (bottom panels) con- 
figurations. Due to the intrinsic noisy nature of the bispec- 
trum, we have binned the values in fc— space, in the same 
way as the flux power of the previous subsection. 

One can see that while at z = 3 the differences are very 
small and usually less than 3-4%, they become much larger 
and of the order of 30-40% at 2 = 4. At z = 5.5 the dif- 
ferences become again smaller and with different wavenum- 
ber dependence. It is possible to interpret this trend in the 
fr amework o f the se cond order perturbation theory as done 
m I Viel et al.l(|2004bh : the overall amplitude and shape of the 
flux bispectrum could not be smooth and strongly depend 
(in a non-trivial way) on the redshift evolution of the coeffi- 
cients that describe the evolution of the mean flux level and 
of the IGM temperature-density. 
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Figure 5. Ratio between the simulated ID flux power spectrum of four different (60, 384) models with /nl = —200, —100, +100, +200 
and /nl = 0, represented by the blue dashed, blue continuous, black continuous and black dashed lines, respectively. Results are shown 
at at z = 3, 4 and 5.5 in the left, middle and right panels, respectively. 
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Figure 6. Ratio between the simulated ID flux bispectrum in the flattened (k, k, —2k) (bottom row) and squeezed configuration (top 
row) (k — k m i„, —k — k m i n , 2k m i n ) with k m i n = 2tt/L (and with L the linear size of the box in km/s) of four different (60, 384) models 
with /nl = —200, —100, +100, +200 and /nl = 0, represented by the blue dashed, blue continuous, black continuous and black dashed 
lines, respectively. Results are shown at z = 3,4, 5.5 in the left, middle and right panel, respectively. 
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Figure 7. Ratio between the simulated flux probability distribution function of two different (20, 256) models with /nl = —200, +200 
and the gaussian model with /nl = 0, represented by the blue and black lines, respectively. The continuous lines are for an infinite 
signal-to-noise ratio, while dashed and dotted are for signal-to-noise ratios of 100 and 250. Results are shown at z = 3, 4, 5 in the left, 
middle and right panels. The continuum errors are modelled as described in the text and are of the order 2% at z = 3 and 10% a t 
2 = 4, 5.5. In the left panel we report as a filled area the statistical error bars of the flux PDF at z = 2.94, as found in l lKim e t al. 2007). 



4 DISCUSSION 

Among the different flux statistics that we have explored, the 
flux PDF seems the most promising in order to detect pri- 
mordial non-Gaussiantity. However, to assess whether such 
information can actually be extracted from the real datasets 
one needs to compare the expected signal with the ampli- 
tude of the known errors. The present statistical uncertain- 
ties in the flux PDF at z < 3 derived from high-resolution 
high signal-to-noise spectra is below 4-5%. This number was 
derived using jack-knife estimators from a suite of high- 
resolution high sign al-to-noise ( > 50 and usually around 
100) QSO spectra bv lKim et all l|2007l ). This is smaller than 
the effect we are seeking and possibly the NG signature is 
degenerate with other effects such as a change in the tem- 
perature evolution of the IGM. 

In Figure [7] we show in a quantitative way the effect of 
the observational errors on the flux PDF, at z — 3, 4 and 5.5. 
We use a r ealistic (observed ) array of signal-to-noise values 
taken from lKim et al.l l|2007l ) at z = 2.94. The signal-to-noise 
depends on the transmitted flux, and the average signal-to- 
noise value for the noise array taken turned out to be ~ 100. 

The various curves in these plots are the same as in 
Fig. [3] for the /nl = ±200 cases only. The dashed refers 
to the realistic errors of|K im et al.l (|2007| ) corresponding to 
an average S/N = 100. We also plot the case of a more 
favorable case with S/N = 250 (dotted curve), while the in- 
finite signal-to-noise error is represented by the continuous 
line. A second source of uncertainty is represented by con- 
tinuum fitting errors that we have modelled in a statistical 
way that produces a ±2%, ±6% and ±10% displacement of 
the continuum level at z — 3 and at z — 4 and z — 5.5, 
respectiv ely. These numbe rs hav e been derived by t he esti- 
mates of iKim et al.l l|2007l ) and iBecker et all (|2007l ) based 
on the analysis of high resolution and high signal-to-noise 
QSO spectra. To account for these errors we have adjusted 



the simulated continuum of the transmitted flux along ev- 
ery line-of-sight by factor 1 ± G, where G is a number drawn 
from Gaussian distributions with width 0.02, 0.06 and 0.1 
at z = 3, 4 and z — 5.5. This should provide a reasonable 
estimate of the continuum fitting errors effects on the flux 
PDF as long as these errors are Gaussian. 

Of course, taking into account realistic signal-to-noise 
values and the continuum fitting errors reduces the signif- 
icance of the NG signal. We find that for a signal-to-noise 
ratio of 100 (250) it is reduced by 40 % (20%) at z = 3, 
in the same way at higher redshifts, where the NG signal is 
higher, we find similar values. The continuum fitting errors 
are somewhat more important and reduce the significance 
of the NG signal on the flux PDF for / NL = ±200 by ~ 40% 
at 2 = 3. However, adding the two sources of errors at the 
same time as shown in Figure [7] decreases the NG signal by 
45 %, 40% and 80% at z = 3, 4, 5.5 for the / NL = ±200 
cases, respectively. 

The statistical errors estimated bv lKim et al.l (|2007t ) are 
represented by the shaded area in the leftmost panel and 
refers to z = 2.94. Ideally, one would like the NG signal to 
be larger than the statistical errors once all the systematic 
errors have been taken into account. At z = 3 we are indeed 
in this case, but only marginally so. We find that the effect 
of including continuum uncertainties at z = 3 has the same 
quantitative effect of dealing with a signal-to-noise ratio of 
250 instead of an infinite one, and when these two errors are 
added together the effect on the flux PDF is of the same 
order of the NG signal for /nl = ±200. At higher redshifts 
the situation becomes slightly better. However, despite the 
reduction of the NG signal, its signature is still large enough 
to be detected, especially at z > 3, and a higher significance 
could be of course reached once all the Lyman-a flux statis- 
tics (PDF, flux power and bispectrum) will be fitted at the 
same time. 

We stress that our quantitative arguments do not in- 
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elude the possible degeneracies on the flux PDF of NG 
with othe r cosmological and a strophysical parameters as ad- 
dressed in lBolton et al.l l|2008l ). It is however intriguing that 
a better fit to the PDF data presented there at F > 0.8 
would require emptier voids and thus negative values of /nl- 



5 CONCLUSIONS AND PERSPECTIVES 

In this work we have explored the possibility of constrain- 
ing primordial non Gaussianity through the statistical prop- 
erties of Lyman-a forest QSO spectra at 2 ^ 3. For this 
purpose, and for the first time, we have performed a suite 
of high resolution non-Gaussian hydrodynamical simula- 
tions. Although recent analyses have provided convincing ev- 
idence that the most stringent constraints to primordial non- 
Gaussianity will be likely provided by t he large scale biasin g 
properties of rare, massive objects (e.g. lSlosar et al.l ((2008)), 
the analysis of the Lyman-a forest has to be regarded as 
complementary since it would probe non-Gaussianity on 
smaller scales and at intermediate epochs between other LSS 
probes and the CMB. 

The main results of this study can be summarized as 
follows: 

(i) the differences between the Gaussian and non- 
Gaussian scenarios are more evident in regions of high flux 
transmissivity associated to low density environments in the 
gas distribution; 

(ii) deviations from the Gaussian case are best seen in 
the high flux tail of the ID flux PDF: differences are of the 
order of 20-30 % and 2 = 3 and increase up to ~ 100 % at 
2 = 5.5; 

(Hi) differences in the void distribution function are 
comparatively smaller, indicating that the PDF is a better 
statistics to spot primordial non-Gaussianity; 

(iv) the ID flux power spectrum is little affected by non- 
Gaussianity, as expected by the analogy with the matter 
power spectrum: the measured differences are of the order 
of a few per cent and increase at higher redshifts; 

(v) the flux bispectrum represents a much more pow- 
erful statistics and potentially could provide strong con- 
straints; 

(vi) the significance of the non-Gaussian signal is highly 
reduced when one accounts for realistic signal-to-noise val- 
ues in the measured flux PDF and continuum fitting errors 
at high redshifts; nevertheless, significant constraints on the 
non-Gaussianity can still be extracted from the analysis of 
the high flux tail of the flux PDF. 

The statistical error bars on the flux power as measured 
using the SDSS data release 3 bv lMcDonaldl (|2006h are usu- 
ally in the range 3-10 % (going from the small scales 0.01 
s/km to the largest 0.001 s/km) in the range 2 = 2 — 4, 
so the NG signal in this case is smaller than the statistical 
error (even though combining all the data points the error 
will be of the order on the power spectrum amplitude will 
become 0.6% and on its slope of ± 0.005). The SDSS data 
release 3 is based on a sample of 3035, increasing the number 
of observed QSO spectra will reduce further the statistical 
error by a factor U Nqso making the NG signature more ev- 
ident, once the degeneracies with all the other cosmological 
and astrophysical parameters will be properly addressed. 

Regarding the flux bispectrum, the present statistical 



error bars at z ~ 2 are of the order 50% IViel et all (|2004bl ). 
as derived from high-resolution spectra, a value that is much 
larger than what is expected from a non-Gaussian signal at 
that redshift, while this value is comparable to what could 
be seen at z ~ 4. Even in this case, in order to study in a 
precise way putative NG signatures in the flux bispectrum, 
more work is needed to address numerical convergence of 
the flux bispectrum and to incorporate the relevant physical 
processes that can affect its shape and amplitude down to 
smaller scales than those probed by the flux power. 

The statistical error bars derived from present data sets 
of QSO spectra at high resolution in the flux PDF function 
are usually below 5% for high transmissivity regions. This 
value is basically determined by the signal-to-noise ratio of 
the spectra and, at least potentially, higher signal-to-noise 
ratios can be achieved and beat down this statistical error. 
This statistics seems promising due to the large number of 
QSO spectra available and to the better understanding of 
systematics. Among the possible systematics the most im- 
portant is the uncertainty due to the continuum fitting er- 
rors which however could probably be significantly reduced 
at high redshifts with a better understanding and removal 
of the QSO continuum. 
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